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Abstract 

Two distinct models for self-similar and self-affine river basins are numerically 
investigated. They yield fractal aggregation patterns following non-trivial 
power laws in experimentally relevant distributions. Previous numerical esti- 
mates on the critical exponents, when existing, are confirmed and superseded. 
A physical motivation for both models in the present framework is also dis- 
cussed. 



PACS numbers:68.70.+W,92.40.Gc,92.40.Fb,64.60.Ht 
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I. INTRODUCTION 



Experimental analyses of river networks |]J have shown clear examples of behavior analo- 
gous to critical phenomena characterized by the absence of a single well-defined length scale 
reflected in a power-law behavior of various quantities. A fundamental question that arises 
from these observations is whether, in analogy with conventional critical phenomena, one 
may fruitfully classify this behavior into universality classes that are characterized by differ- 
ent sets of exponents. A related point is whether there exist scaling relationships between 
the various exponents of a given universality class. Another vital issue is the elucidation of 
simple models amenable to analysis that nevertheless capture some of the features of fractal 
fluvial patterns. 

Within this framework a theoretical description of the system needs to address two 
basic issues [@,||,@,^@J^^| • First, a careful characterization of the topological properties 
of the networks is essential for understanding the basic transport mechanism in the basin. 
References 0H are recent attempts in this direction. Optimization principles have been 
exploited both numerically 0,|7j and analytically || to explain the tendency of natural 
drainage networks to evolve toward an optimal stable topology. General scaling arguments 
can be found in || . 

Second, a study of the dynamical evolution of the landscape (on geological time scales) 
as a result of interaction with external agents (rain, wind etc) would be desirable HH. 

The present work will address only the first point. We will argue that, despite much 
progress in the past few years, the problem is not yet fully understood and deserves further 
analysis. To this aim we will discuss, based on physical arguments, two toy models of 
river networks. The first model leads to a self-similar river basin, and is relevant when 
the erosional properties of the surface soil are strongly heterogeneous. The second model 
considers the homogeneous basin case and results in a self-affine river network. 

Although the overwhelming majority of the observational data are consistent with a self- 
affine description (i.e. networks display a privileged direction), the marked self-similarity of 
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the basins with their own sub-basins suggests a crossover from a self-similar character above 
some length scale. This is one of the reasons for considering models with both characters. It 
should also be emphasized that although some of the features of the models presented here 
were previously discussed in the literature (see below), we believe that both the physical 
motivations and the analysis carried out here are essentially new. 

The plan of the paper is as follows. In the next section few definitions and scaling 
relations will be recalled. In Sec. Ill results for self-similar are presented. Sec. IV is 
dedicated to the self-affine counterpart. In Sec. V few relevant experimental results will be 
briefly reminded for the sake of completeness. Sec. VI will summarize our findings along 
with some future perspectives. 



II. DEFINITIONS AND SCALING LAWS 



We define a river network as a spanning (loopless) tree on a lattice of linear size L |L0 
Each site has exactly one output bond to one of its neighbors and no restriction on the 
number of input bonds (three at most on a square lattice). In a river basin, the area at any 
site is defined as the number of sites upstream of the site connected by the network. From 
the computational point of view, it can also be regarded as a measure of the flow rate if a 
unit weight is assigned to each source thus simulating a unit constant precipitation. 

The equation for s i5 the area at a given site i, is 

j£nn(i) 

where Wij is 1 if i collects water from its nearest-neighbor (nn) site j and otherwise. 

It is experimentally observed |IJ and theoretical explained [^,|J that in river basins the 
probability density p(s, L) of a site having area s in a system of size L, has the scaling form 

p(s,L)=s-F(^) (2) 

where F(x) is a scaling function which takes into account finite size effects and is the finite 
size exponent. 



Similarly the distribution of upstream lengths has been also predicted || and confirmed 



by field observation [11] to display the universal form: 



<i,L) = ry(^). (3) 

where f(x) is the analogue of F(x) and d\ coincides with the stream fractal dimension. 
The upstream length is defined as follows. At a given site the areas (see eqn. ([[])) of the 
nearest-neighbours are checked. The site with the largest value leads to the outlet. The site 
with the next-largest value is defined to be an upstream site - it indicates the longest path 
towards the source. If two (or more) equal areas are encountered, one is randomly selected. 
Alternatively, a burning algorithm Q could be also employed. 

In natural basins, the drainage area s and the stream length I are related by Hack's law 







s ~ l x ' h (4) 

The sub-basin from any site defined as all the upstream sites connected to it is character- 
ized by typical longitudinal and transverse lengths £y and For self-affme river networks 
one defines the Hurst (or wandering) exponent as £j_ ~ Cjf with H < 1. Note that for 
self-similar river networks (in which each rivulet originating from any site and proceeding 
to the global outlet is a fractal characterized by the same fractal dimension di), the Hurst 
exponent H — 1. 

As one might expect the exponents are not independent. For self-similar networks (di > 
l,H = 1) di determines all the other exponents ||, 

= 2 , h = d l /2 , r = 2-2/d l , 7 = 2M (5) 

For self-affine networks (di = 1,H < 1) it is H which defines all the other exponents p] 

4>=1 + H , h = l/(l + H) , r = ^§- , j=l + H (6) 

1 + ti 

Two features of the above relations are worth mentioning. First there is experimental 
evidence in the observed data that H < 1 and di > 1. This apparent contradiction might 



be explained with the crossover between the two regimes occurring at some length scale, 
as mentioned in the introduction. Secondly it turns out from (|5]) and @ that identical 
values of the exponents are obtained from both cases if d\ = 2/(1 + H). This means that 
knowledge of the exponents other than d\ and H cannot discriminate the self-similar or 
self-affine character of the basin. In this respect a direct measure of di and H appears to be 
crucial for its characterization. 



III. SELF-SIMILAR RIVER NETWORK MODEL 

We first discuss the model of self-similar river networks. Consider a network which is a 
square lattice of size L x L where the links of the rivers are identified with the bonds of the 
lattice. Periodic boundary conditions are assumed in the left-right direction. The bottom 
side of the square is defined to be the (fixed) outlet which collects the water that is flowing 
out. Independent random numbers in the range (0,1) are assigned to the different bonds 
representing the erodability P[s, of the surface soil of the bond i. 

The physical situation we have in mind leads to river network formation based on an 
invasion percolation like mechanism fl3| . The weakest erodable link is selected and assumed 



to be a part of the network. The second-ranking weakest link is then selected and so on. 
The process is iterated in the ensemble of the remaining links until all sites are connected, 
i.e. they all have a route to the outlet. Loops are excluded since once a preferred route 
is selected, alternative routes formed due to the presence of a loop would be energetically 
unfavourable. Operationally, one thus obtains the network by incorporating the regions in 
order of increasing strength so that no loops are formed, yet all sites on the lattice are con- 
nected to the outlet sites. A variant of the above procedure leading to the same structure, 
consists of starting from the links connected to the outlets, selecting the weakest one and 
proceeding invasively (i.e. always choosing the new weakest link) in the new ensemble of 
the interfacial links. This model, which was originally introduced by Stark |2| was subse- 
quently rediscovered by Manna and Subramanian B. This is a model of headward growth 
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of streams away from a rift, the weakest bond corresponding to the point most susceptible 
to bank failure. The motivations which led the above authors to the model were however 
completely different from ours. One might suspect that a variant of the above model having 
(statistically) spherical geometry would lead to a different different universality class. We 
checked that this is not the case by starting with a central outlet and proceeding as above 
until the whole domain is drained. 

A typical river obtained by our procedure is shown in Figure |l|. We have carried out 
detailed studies of the scaling properties of the networks. Our numerical simulations involved 
sizes up to L = 192 with a typical number of different configurations of the order of 500. 
A summary of our results is presented in the Table along with the results of observational 
data |l|]. In order to get a precise estimate of the r exponent we used two different methods. 
The first consists in plotting the local slope (which is trivially related to r) of the cumulate 
area distribution (see Fig.0) 

P(s,L)= f ds'p(s',L) (7) 
Jo 

as a function of s. This is reported in Fig.|3]. An average over all these values yields 
t = 1.406 ± 0.021. On the other hand we performed a finite size analysis of the exponent 
extracted at the various sizes (see Fig.f|). An extrapolation then yields r = 1.404 ± 0.001. 
The value reported in the Table is then the arithmetic average of these two values. The 
exponent <fi can be determined by plotting the universal function as defined in (0). This is 
also shown in Fig.|5]. In a similar we computed the exponent 7 as defined in (|3|) obtaining 
7 = 1.612 ± 0.049. In Fig. || the universal function f(x) is computed yielding a good 
collapse for d\ = 1.22. The value for di can be confirmed by an independent computation of 
the typical length 

<r l q > 

= L dl (8) 

where the averages are referred to the probability distribution density ir(l, L). Using various 
values of q we found stable values allowing an estimate as di = 1.22 ± 0.05. 



The scaling predictions are found to hold very well. Our exponents agree with those 
numerically obtained recently in Ref. || Specifically they reported r = 1.392 ± 0.010 and 
7 = 1.628 ± 0.05, upon using sizes up to L = 1024 but with a less precise data analysis. 

This model may be alternately viewed as an optimization problem that selects the span- 
ning tree that minimizes J2i Pi where Pi denotes the erodability of the i — th bond and the 
sum over i runs over all the bonds of the tree. [f|[n|] Our model then becomes a special case 
of the recently introduced optimal channel networks fijjTj that are constrained to satisfy a 
global minimization of the energy expenditure. |16 



We note that the resulting network is a union of the invasion percolation paths from each 



of the sites to the outlet |[17|| . We also note that when two paths intersect, they coincide 
the rest of the way to the outlet. This provides a natural mechanism for aggregation in 
our model. Also each of the individual paths corresponds to the optimal path in a strongly 
disordered medium that was recently shown |18|] to be characterized by di = 1.21 ± 0.02. 
Note that this value is somewhat different from the value 1.13 cited by Stark [Q]. That 
value corresponds to the fractal dimension of the geometrical shortest path on a percolation 
cluster and is different from the fractal dimension of the loopless strands which correspond 
to the energetically shortest path in the sense explained above. 

The Hurst exponent H obtained |19| from a box-counting dimension d\ would be 2 — d\ 



and thus in the range 0.77-0.81. With this value of H and using the scaling relations in the 
text, we obtain r = 1.44 ± 0.04, h = 0.56 ± 0.01, and 7 = 1.79 ± 0.02 in perfect agreement 
with the observational data. Such a crossover might occur at some intermediate length scale 
beyond which the rivers are no longer self-similar. 

We note that the homogeneous analog of our optimal channel network (Pj=const) is the 
random spanning tree problem, in which all trees occur with equal probability, d\ = 5/4 
and the other exponents follow from our scaling relationships [|10j. Our model seems to be 
in a different universality class presumably due to the distinct weights associated with the 
different trees. The possibility that our result of d\ » 1.21 may crossover to 5/4 for much 
larger sizes cannot of course be ruled out from our data [20 



IV. SELF-AFFINE RIVER NETWORK MODEL 



We now turn to the second model leading to self-affine (H < 1) river networks. This 
corresponds to the homogeneous version of the first model where the Pj's are all equal. Now 
the interfacial links all have an equal probability of being invaded with the usual constraint 
that loops are not formed. This procedure is akin to the well-known Eden growth problem. 
A loopless cluster generated in a two-dimensional Eden growth process on a square lattice 
with a central seed site is shown in Figure [7]. This model is not new. Meakin |HJ numerically 
studied the same model although with different aims. By mapping this problem onto a 1 + 1 
Kardar-Parisi- Zhang equation, Krug and Meakin [^] realized that the value r = 1.40 should 



be exact. As argued in Refs p3fl , the individual rivers are no longer self-similar but self-affine 
with a Hurst exponent of 2/3. The Hurst exponent is larger than the random walk value of 
| because the Eden growth process generates strands that compete with each other p3[ and 



effectively mimic a quenched disordered environment || . The exponents predicted from this 
value of H are shown in the Table. Figure |] shows the log-log plot of P(s, L) vs. s. The data 
points analyzed along the same lines as before, are consistent with the exponent r of 1.40 
which agrees with the scaling prediction (see Table). It is remarkable that the self-similar 
and self-affine river network models yield very close predictions (indeed di ~ 1.21 whereas 
2/(1 + H) =6/5) for all of the scaling exponents even though the underlying mechanisms 
are very distinct. We believe this point deserves further investigation. 

V. OBSERVATIONAL DATA 

In this section we shall review some experimental known results to provide evidence that 
the two above toy models, albeit rather primitive, are in fact relevant in the interpretation 



of observational data. For a more exhaustive analysis see Ref. ||TT . 

The network associated with a given natural terrain pertaining to a river basin can be 
experimentally analyzed by using the so-called Digital Elevation Map (DEM) technique 
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(see e.g. and references therein). In its digitized form, the elevation map of a terrain 
allows the determination of the soil height of areas (pixels) of order 10 -2 Km 2 . Thus a 
fluvial basin is represented in a objective manner over few (typically 3-4) log scales of linear 
size. A flowrate unit is associated with each pixel and the flow contributing to any pixel 
follows the steepest descent path through drainage directions. The resulting network (thus 
defined by the drainage directions) is therefore a two-dimensional representation of the 
three-dimensional landscape. 

In Fig. P a representative network obtained in this fashion is shown. This particular 
network has the exponent r = 1.40 which is in very good agreement with both the self-similar 
and the self-affine models. 

It is important to stress that the actual values of the critical exponents do vary from 
basin to basin. However, within each single river network, all the exponents closely satisfy 
the scaling relationships that we have derived. The exponent values for the network of Fig. 
|9] is in perfect agreement with our self-affine (Eden-like) model, although other networks 
often have slightly higher values of r [IT|. 



The measured values for d\ and H ranges between 1.02 and 1.07 and between 0.75 to 
1.00 respectively ]TT |. The above models then lie at the outer edge of the observed data in 



both cases. We believe that although these two models do represent, at a very schematic 
level, two physical mechanisms occurring, at some length scale, in real rivers, our numerical 
results show that this is not sufficient and that other effect possibly combine to yield different 
numerical exponents. 

As mentioned earlier, identical statistics can be obtained by using either self-similar 
(di > 1, H = 1) or self-affine (di = 1, H > 1) networks, provided that di = In this 

respect a direct measure of the anisotropy of real rivers appears to be crucial. 

Finally we remark that although the observational data favour self-affine rivers, the 
(statistical) similarity between each basin and its sub-basins, suggests that river networks 
are indeed self-similar, although their actual shape is anisotropic (as expected for a self-affine 
network) . 
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VI. SUMMARY 



We have studied two lattice models of spanning trees that lead naturally to fractal river 
networks. We have chosen a self-similar and a self-affine case on the basis that in nature 
there are features belonging to both characters, although it is commonly accepted that 
natural rivers are self-affine. The two models have very different physical motivations. The 
self-similar model is directly driven by disorder. Our choice of heterogeneities is spatially 
uncorrelated, and any degree of spatial correlation would decrease the small scale tortuosity 
of disorder-dominated paths. The resulting spanning tree-like structures exhibit a general 
similarity to river basins in overall appearance and very consistent scaling statistics. Our 
numerical results are consistent with previous estimates of the same model introduced on 



different physical grounds and corrects some misunderstandings present in the literature 
0. More generally, our results show that fractal structures arise from the minimization of 
a disorder-dominated total energy functional and reinforce earlier suggestions [f|] on the 
connections of optimality with fractal growth. On the other hand the self-affine model has 
no disorder in the definition but it has buried in it a competition mechanism which effectively 
yields quenched disorder. Again our results are in accord with previous investigations but 
have broader consequences in our framework of analyzing the effects of heterogeneities on 
the networks. 



Two recent investigations are pertinent to our work. In Ref. |24[ the interplay between 
quenched disorder and non-linearity in the landscape evolution was shown to be relevant 
for the interpretation of the real rivers field data. On the other hand, a even more recent 
renormalization group analysis of a continuum equation p5 suggests that the two disorder- 
dominated networks studied here, the unweighted spanning-trees studied in Ref. and the 



aforementioned landscape model of Ref. ||24|| , all belong to different universality classes albeit 
with very close exponents. We believe that this scenario is intriguing and deserves further 
attention. 
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Table 





SELF-SIMILAR 




SELF-AFFINE 








Scaling Predictions 


Measured 


Scaling Predictions Measured 


River Basins 




(with di = 1.21 ±0.02) 




(with H = §) 






di 


1.21 ± 0.02 


1.22 ± 0.04 


1 


1 


1.1 ± 0.02 


H 


1 


1 


2 

3 




0.75 - 0.80 


T 


1.395 ± 0.01 


1.38 ± 0.03 


7 
5 


1.40±0.02 


1.43 ± 0.02 


h 


0.605 ± 0.01 


0.62 ± 0.02 


3 
5 




0.57-0.60 


7 


1.65 ± 0.03 


1.60 ± 0.05 


5 
3 




1.8 - 1.9 


dp 


1.21 ± 0.02 


1.21 ± 0.02 


1 


1 









The exponents predicted by the scaling arguments, measured in our simulations and for 
river basins, dp is the fractal dimension of the river basin boundary. Note the inconsistency 
in the observational data - d\ is greater than 1 suggesting a self-similar network, whereas 
H < 1 indicating a self-affine structure. 
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FIGURES 

FIG. 1. Typical self-similar river network on a 128 x 128 lattice obtained by our optimization 
procedure. Only the largest river is shown. Periodic boundary condition are used only in the 
direction transverse to the dominant flow. The size of the circles is a measure of the value of s. 

FIG. 2. A log-log plot of the cumulate area distribution P(s, L) vs. s, for lengths ranging from 
L = 16 to L = 96. Two slopes corresponding to the Scheidegger (dotted line) and to the Mean 
Field model (full line) [l^J] are also shown as a guide for the eye. 

FIG. 3. Effective exponent r as computed from the local slope. 

FIG. 4. Plot of the exponent r as a function of 1/L. The extrapolated value is r = 1.404±0.001 
is indicated by a star. 

FIG. 5. Collapse of the p(s,L) curves for L=8 (solid diamonds), 16 (boxes), 32 (solid 
hexagons), 64 (stars), and 96 ( closed circles) according to eq.2, with r = 1.40 and <j) = 2.0. 

FIG. 6. Plot of the universal function f(x) obtained from eqn (|| for the same values as in 
Fig.5) with 7 = 1.61 and di = 1.22. 

FIG. 7. Typical network of self-affine rivers flowing to an outlet in the center. The network 
has been obtained by generating an Eden growth process from a central seed and stopping the 
growth when the maximal horizontal distance reached is equal to 64 lattice constants. The size of 
the circles is a measure of s. Only sites with s > 100 are shown. 

FIG. 8. A log-log plot of the cumulate area distribution P{s,L) vs. s for the self-affine model. 
The system sizes are L=32, 64, and 256 as indicated in the figure. The number of samples for the 
three cases are 1000, 500, and 200 respectively. The solid line corresponds to the exponent r=1.40. 

FIG. 9. The transportation network of the Johns river drainage basin in Kentucky, USA. Its 
extension is 984Km 2 . The measured exponents are: r = 1.40, if = 0.67. 
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